Nonequilibrium Phase Transition in Non-Local and Nonlinear Diffusion Model 
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We present the results of analytical and numerical studies of a one- dimensional nonlocal and non- 
linear diffusion equation describing non-equilibrium processes ranging from aggregation phenomena 
to cooperation of individuals. We study a dynamical phase transition that is obtained on tuning 
the initial conditions and demonstrate universality and characterize the critical behavior. The crit- 
ical state is shown to be reached in a self-organized manner on dynamically evolving the diffusion 
equation subjected to a mirror symmmetry transformation. 



PACS numbers; 05.20.-y 02.50.Le, 05.40.-a, 05.45. -f-a 



The study of equilibrium phase transitions is now a 
mature field ||l|,^ and, increasingly, attention is being 
paid to outstanding problems in nonequilibrium statis- 
tical physics H]. Such problems are often challenging 
because of inherent non-linearities. There are many 
examples of non-equilibrium phenomena which are in- 
trinsically non-local such as the growth of thin films in 
the presence of shadowing |^ and the sculpting of the 
drainage basin of river networks due to erosional pro- 
cesses 1^. A striking development in the field of non- 
equilibrium statistical physics is the development of the 
paradigm of self-organized criticality entailing the 
competition between two dynamical processes leading to 
a critical state without any fine tuning of parameters. 

The principal theme of this letter is the study of a 
simple one dimensional nonlinear and nonlocal diffusion 
equation to elucidate the nature of a non-equilibrium 
phase transition. The equation is essentially one de- 
scribing biased diffusion Q with the magnitude of the 
bias determined by the instantaneous configuration of 
the random walkers undergoing diffusion. Physically, 
the equation describes aggregation, population dynam- 
ics, and represents a simple model for cooperative behav- 
ior in game theory Strikingly, the behavior changes 
from that of a conventional critical point (which requires 
tuning) to that of self-organized criticality Q on con- 
sidering the time evolution of a transformed equation 
obtained by a mirror symmetry transformation with x 
replaced by — x. 

Our basic equation for P{x,t), the probability that the 
diffusing particle is at position x at time t, is 



dP _ dP 1 d'^P 
'dt ~ ^''^^'~d^ ^ 1~d^ 



(1) 



where both the nonlinearity as well as the nonlocality are 
introduced in the bias velocity v defined by 



v{t) = / dxP{x,t) ~ 1/2 . 



(2) 



On setting v equal to zero in Eq. (|^), one recovers the 
standard unbiased diffusion equation whereas 



one obtains simple biased diffusion, when v is a constant. 
In our equation, w is a measure of the imbalance between 
the population of walkers in the right and left and the 
drift bias promotes further aggregation. 

Eq. (|l|) describes the temporal evolution of the distri- 
bution function P{x, t) and leads to one of two outcomes 
in the large time limit. Depending on the initial distribu- 
tion, one ends up with the bias to the right or to the left 
winning so that P becomes 1 either at +oo or at — oo. 
Our focus is on the non-equilibrium transition between 
these limiting behaviors. Note that there is a set of ini- 
tial conditions (of measure zero) that correspond to the 
dynamical critical point - we will demonstrate that the 
critical behavior is universal. 

Let us define a new variable w{t) ~ drv^r), and 
introduce y{t) ^ x — w(t) so that Eq. (|l|) is cast in the 
form of a standard diffusion equation. 



dP _ 1 d^P 
'dt ~ 1'df 



(3) 



For simplicity, let us first consider an initial Gaussian dis- 
tribution of P(x, t) centered around and with variance 
(Tq. The solution of Eq. (^ is then given by 



P{x,t) = iVexp <^ - 



[x~xa- ■w{t)f 



(4) 



where N = 1 / ^/^^(T+2(J^ is the normalization constant. 

With (To 0, P{x,t) is also the fundamental solu- 
tion jl^, which will be used to obtain the distribution 
at time t, starting from more general initial conditions. 
Expression (||) is only a formal solution because w{t) is 
itself a function of P{x, t). The transition is between two 
phases corresponding to aggregation on the right or on 
the left and therefore one would expect that the criti- 
cal point would correspond to a situation with no bias 
(i.e. Xq = 0). In this case, the distribution is symmet- 
ric with respect to the origin at all times, and there is 
nothing to choose between left and right, thus an unbi- 
ased behaviour ensues. In order to probe the behaviour 
in the vicinity of this critical point, one could start with 



1 



an initial distribution with a tiny value of Xq (small bias) 
and watch how the system evolves. Combining Eqs. (|l|) 
and (0) and noting that w — v, we find 



w — —bjTli 

2 



2al 



(5) 



At criticality {xq — 0), the velocity is zero at all times 
and w vanishes too. In the critical regime, one expects 
that w{t) ^ I, for all t less than a crossover time Tr, so 
that one may linearize Eq. to find that 



Xo 



V^t + 2^2 

This equation can be easily solved and yields 



(6) 



(7) 



One may define a characteristic transient time, Tr, spent 
in the critical region during which the linearization ap- 
proximation holds. From ([7|), one finds that diverges 
as Xo ^ as 



wit) = xo SV^o~V^.)/V^ _ 1 



In'lxol 



(8) 



In the critical region, the characteristic length-scale ^ is 
expected to follow the diffusion law, and therefore 



ln|xo| 



(9) 



as we will verify numerically later on. 

Alternatively, in the initial condition, one may in- 
troduce a bias by fixing xq — and instead letting 
P(oo, 0) = 4>o as an effective asymmetric boundary con- 
dition. With the same procedure used to derive Eq. (||) 
we have 



Erf' 



(10) 



= corresponds to the critical point, and by lineariz- 
ing Eq. ( p^ ) for small 4>o, it is straightforward to show 
that one again obtains Tr ~ ln^($o), with ^ ~ \/T^- 

We have verified that the same critical behaviour holds 
for other families of initial conditions. For example, when 

P{x, 0) = (1/2 - e) 6{x - 1) + (1/2 + e) S{x + 1), (11) 

we find that 



Tr 



In' lei 



(12) 



and ^ ^ |ln|e||. Note that e is now a measure of the 
deviation from the critical point and s = correponds to 
the unbiased, zero drift situation. 

We now turn to the results of numerical experiments on 
a 1-dimensional lattice which are useful for validating our 
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FIG. 1. Behaviour of transient time Tr (at which v{t) first 
becomes equal to 1/2), for the initial distribution (|ll]) for dif- 
ferent e values. In the plot, the square logarithmic divergence 
corresponds to linear behaviour. 



analytic predictions and for probing the nonlinear regime. 
The discrete version of Eq. d), used in our simulations, 
reads 



P.(t+1) = 



lit) 



1 - $(i) 



(13) 



with (f>(t) = J2x>o P^i^)- velocity is thus given by 
v{t) = ^{t) — 1/2. This equation was proposed by Nowak 
and Sigmund |lll| as a simplified model of the evolution 
of indirect reciprocity by image scoring. Indirect reci- 
procity is determined by reputation and status and is 
characterized by each individual having an image score. 
A potential donor and recipient of an altruistic act have 
an interaction in which the donor helps the recipient pro- 
vided the recipient's image score is positive. Such an al- 
truistic act increases the image score of the donor by 1 
(the selfish act would have decreased it by 1) and the 
image score of the recipients is unchanged. Eq. (^) is 
the governing equation for the time evolution of P^, the 
fraction of players with image score x. The two phases 
that we have considered correspond to cooperation and 
defection and our finding is that not much time is spent 
agonizing over which phase to select even in the vicinity 
of the critical point. Indeed, the time scale to decide on 
one of the two different phases only diverges logarithmi- 
cally as one approaches the critical point. 

Fig. shows the divergence of as e — > for the 
initial distributions given in Eq. ( pi] ) and for a lattice 
of 6x10^ sites. Numerical results are in excellent accord 
with our theoretical predictions (see Eq. (|l2|)) 
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FIG. 2. Temporal behaviour of the bias velocity v{t) for 
various values of the parameter $o around the critical value 
$c. The upper curves correspond to values of $o approaching 
$c from above. The lower curves are obtained when $o tends 
to from below. 



FIG. 3. Dependence of the transient time Tr on ($c — "I'd). 
Tr is defined as the time needed for v{t) to approach the 
value 1/2 when <I>o > and $o — 1/2 for <I>o < <l>c. The 
bullets show the behavior when $o < while the triangles 
correspond to $o > *&c. 



Let us consider now, an interesting generalization 
of the initial conditions that we studied analytically: 
P(oo, 0) = $0 and 0) = 1 - $o and P = at 

all other locations, initially. The critical value of $0 in- 
creases monotonically to a non-zero value ^dz) for pos- 
itive z with ^c{z 00) = 1/2. The smallest value of 
$c(^) on a discrete lattice occurs when z = 1, which 
is the case we focus on. $c = ^c{z = 1) is found to 
be 0.261970531164..., a result that was noted earlier by 
Nowak and Sigmund |pT|] . 

Fig. H shows the behaviour of the bias velocity as a 
function of time as $0 approaches $c from above and 
from below. Tr can be identified as the time after which 
v{t) becomes equal to its asymptotic values of either 1/2 
and $0 — 1/2 and its scaling behaviour is shown in Fig. ||. 
The average location of the random walkers (excluding 
the number fixed at a; = cxd) behaves with time (in the 
vicinity of the fixed point) as 

{x{t)) ^ .y< + 2^ . (14) 

In order to derive this result, we note that the velocity 
increases very slowly in the linear critical regime and can 
be approximately considered constant (this is consistent 
with the assumption that, for time i <C TV, w is much 
smaller than 1). The derivative of the velocity given by 
Eq. (^ is then vanishing, and using Eq. (|l|) to eliminate 
P we obtain, after integrating over x, the expression 

AvP{0,t)-d^P{0,t) = . (15) 



Using the formal solution (Q), we obtain 

w{t) - \lt + 2al 

The result (|l^) then follows on noting that the average 
position of the walkers is given by w(t) + xq. 

This critical regime behaviour of {x{t)) crosses over 
to a linear temporal behavior when the bias reaches a 
sufhcient strength. There is indeed a sharp onset of the 
linear behaviour at a value of {x) , which one may identify 
with ^. The scaling behaviour of ^ is shown in Fig. ^. 

We now turn to a simple mechanism for obtaining self- 
organized critical behavior in our model. In Eq. (|l|) the 
transformation x — > —x is equivalent to a change of the 
sign of the bias velocity. In this situation, the system 
spontaneously organizes in such a way that the aggre- 
gation of walkers is disfavoured. As a consequence, the 
asymptotic distribution becomes symmetric (character- 
ized by u = 0) and this corresponds to a self-tuning to 
the critical state, a behaviour typical of self-organized 
criticality 

In this asymptotic regime, the scaling {x{t)) ~ ^/t de- 
rived in (|lj) still ought to hold and is confirmed by sim- 
ulations performed on the discretized diffusion equation 
(see Fig. I). 

In this letter, we have introduced and studied a diffu- 
sion equation with nonlinear and nonlocal features. In 
this model, spontaneous fluctuations in the population 
of walkers are able to drive the entire population to- 
wards one of the two boundaries located at (x = ±00). 
This mechanism is of interest as the basis for the devel- 
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FIG. 4. Plot of the characteristic length ^, as a function of 
the deviation from criticality "I>c — $o- This numerical result 
confirms that, in the critical regime, ^ diverges logarithmically 
as the critical point is approached. 

opment of more realistic models of self-aggregation and 
self-organization in cooperative states of populations of 
interacting individuals. A mirror symmetry transforma- 
tion applied to the equation reveals a dynamical evolu- 
tion corresponding to generic behavior associated with 
self-organized criticality. 
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